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Abstract. The work analyzes a one-dimensional viscoelastic model of blood vessel growth 
under nonlinear friction with surroundings, and provides numerical simulations for various 
growing cases. For the nonlinear differential equations, two sufficient conditions are proven 
to guarantee the global existence of biologically meaningful solutions. Examples with break- 
down solutions are captured by numerical approximations. Numerical simulations demon- 
strate this model can reproduce angiogenesis experiments under various biological conditions 
including blood vessel extension without proliferation and blood vessel regression. 

1. Introduction 

Angiogenesis, the growth of new blood vessel from pre-existing vasculature, is crucial 
to many physiological and pathological processes, such as embryonic development, tumor 
growth, and wound healing. Typically these new vessels are very thin (thus also called 
capillaries) and are lined up by an array of tightly adhered endothelial cells. Endothelial 
cells migrate up the gradient of chemotactic cues such as vascular endothelial growth factor 
(VEGF) released by tumor cells. Along a capillary, the tip cell generates a protrusion force 
to drag the whole capillary forward, and at the same time, stalk cells after the tip proliferate 
and generate new masses to sustain the capillary extension. During the capillary extension, 
endothelial cells have to overcome the adhesion (or friction) with surroundings, such as other 
types of cells and extracellular matrix. A detailed description of angiogenesis biology can be 
found in [I]. 

Traditional mathematical models on the capillary extension focus on the chemotactic 
migration of endothelial cells which is typically modeled as a reaction-convection-diffusion 
equation (see reviews in [I] and [2]). Such models lack the investigation of relationships be- 
tween the growth of endothelial cells, mechanical extension of the capillary, and interactions 
between the capillary and surroundings. These topics are central not only to angiogenesis 
but also to the biological growth of all soft tissues. However, "the mathematical modeling of 
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biological growth is currently in an adolescent stage" ([3]), because many constitutive rela- 
tions are still unknown. Therefore, it is a big challenge to model the relationships mentioned 
above. 

Among efforts to meet this challenge are one-dimensional models of epithelial cell mono- 
layers migration ([I],[S]) in wound healing, and a one-dimensional viscoelastic angiogenesis 
model proposed in [6]. Common features of these models include: (i) the blood vessel cap- 
illary or epithelial cell mono-layer is regarded as viscoelastic material connected by tight 
cell-cell junctions; (ii) the cell growth is treated as an isotropic pressure; (iii) the nonlinear 
friction exists between cells and surroundings. In this paper, we will first analyze the angio- 
genesis model in [6], and then utilize a numerical scheme to find approximate solutions to 
this nonlinear problem. 

Assume, initially, a blood vessel capillary occupies the entire interval [0, 1], and every point 
of the capillary is denoted by the Lagrangian coordinate, x, where x = denotes the root 
which is fixed in space, and x = 1 represents the tip which is free to move. The deformation 
of the point x at time t is denoted as u(x,t). The non-dimensionalized model is 

du da 



du 
dx 



(1-1, 



where (3 is the friction between endothelial cells and surroundings, a is the stress tensor 
defined by 

du du , „. . . 

ji is the viscosity of cells, u = ^ is the material derivative, and f(x, t) is the endothelial cell 
mass density. The initial and boundary conditions are 

u(x,t = 0) = 0, (1.3) 
u(x = Q,t) = 0, a(x = 1, t) = g(t). (1.4) 

where the function g(t) is the protrusion force exerted by the capillary tip cell. 

Note that (l+u x ) is the deformation rate of endothelial cells in the capillary. The equation 
( II. ip has the absolute sign of (l+u x ), so it may admit solutions of negative deformation rate. 
However, a negative deformation rate in one-dimensional space is not biological: it means 
one cell simply passes through another cell. Therefore, we only seek biological solutions as 
defined below. 



Definition 1.1. Biological solution. If (1 + u x ) > 0, then u is a biological solution. 
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To penalize the non-biological solutions, we remove the absolute sign of (l+u x ) in Eq. (11. ip . 
For simplicity, we assume (3 = 1 and \i = 0. The justification to set \i — is provided in 
§ 13.21 and the Appendix. Therefore, the problem we will consider in the main body of this 
paper becomes 



Remark 1.1. If(l + u x ) happens to be negative, then the first equation of (11.51) would behave 
like a backward heat equation, which is notorious for the ill-posedness of the problem. 

Remark 1.2. We will see that the endothelial cell density f(x, t) is dominant in the existence 
of biological solutions. Because its presence is similar to the pressure term in the stress tensor 
of fluid mechanics, we will also call it the pressure. 

The theoretical result of this paper is stated in Theorem 2.1 of §[2} which shows that under 
certain assumptions on / and g, the biological solutions exist globally. The key technique 
to prove Theorem 2.1 is to find some nontrivial transformations to simplify the problem for 
which the maximum principle can be applied. Because of the nonlinearity, it is hard to find 
the explicit solutions. Therefore, in §[3] we develop a numerical scheme to solve the system 
(I1.5P and use it to examine the conditions in Theorem 2.1, including the cases when global 
biological solutions exist and solutions break down in finite time. In §HJ three examples of 
angiogenesis under different biological conditions are presented using numerical simulations. 
Finally, we summarize our analysis in §[5j An appendix gives some descriptions of the model 
and the discussion of its parameters. 

2. Analysis of solution existence and asymptotic behavior 

Before solving the problem, we investigate the compatibility conditions of the initial and 
boundary conditions. It follows from the initial condition in (11.51) that u xx = at (x, t) = 
(0, 0). On the other hand, the boundary condition at x = in (11.51) yields that u t (0, 0) = 0. 
Using the differential equation in (jl.5p . we see that / must satisfy f x (0, 0) = 0. 



u(0,t) = 0, —(l,t) = f(l,t) + g(t)-l, 
k u(a;,0) = 0. 




(1.5) 



Theorem 2.1. Suppose f £ C 2 , g £ C 1 



f x (0,0)=0, f>0,g>0, and 



mi(f(l,t)+g(t))>0. 



(2.1) 
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(a) Suppose f x (0,t) = and f xx < 0, then there exists a global biological solution for 
the problem U.5\) . If, moreover, f(x,t) = f(x) and g{t) = g (nonnegative constant) 
when t > T for some T > 0, then 

\u{x,t)-^{x)-&{x)\<Ce- at (2.2) 

for some a > 0, where 

&{x) = [ (f(s) - l)ds, <g{x) = gx. (2.3) 



(b) Suppose f(x,t) = f(x), then there exists a global biological solution for ( fi.5j) provided 
that 

max fix) — min fix) < 1. (2.4) 

xe[o,i] xe[o,i] 

//, in addition, g(t) = g (constant), then the solution of the problem U.5\) also satisfies 



Remark 2.1. In part (a) of the Theorem \2.1l the assumptions f x (0) = and f xx < 
indicate that f is a decreasing function. One such an example occurs when the endothelial 
cells are dying from the tip, so the capillary regresses. In part (b), the assumptions that f is 
steady in time and its variation is less than unity imply the cell density along the capillary 
is fixed in time and is quite uniform in space. 



Proof of Theorem Set F(x, t) = f(x, t) - 1. 
(Part a). Set 

cu(x, t) = x + u(x, t). 
Then it follows from (jl.5p that u satisfies the equation 

U x U t = LO xx - F x (x,t). 

and 

u(0,t)=0, -£(l,t) = f(l,t) + g(t), 
u(x, 0) = x. 

Therefore, 

^ra F x {x, t) 



U! x (jj x 
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Because F x (0,t) = and u(0,t) = 0, we have u xx (0,t) = 0. If we extend the functions as 
follows 

w(x,t), if < x < 1, f F(x,t), if <x < 1, 

£=<J F(x,t) = < 

- w(-x,t), if - 1 < x < 0, I F(-x,t), if - 1 < x < 0, 



then G C 2 ([0, 1]) and satisfies 



ut = ^ - ^kfe^l. (2.5) 



Differentiating equation (12.51) with respect to x, one has 



{&x)xx (wxi) 2 F xx (x,t) t F x {x,t)u) a 
Let (f(t,x) = Cj x . Then satisfies 



(^x)t — — z jz x 2 ~ 1" 



V^xx V^a: , Fx'-Px F xx 
<Pt= 2 + 2 • 2 - 6 



If, in addition, F xx < 0, then 



. <£xx {fx - Fx) 

ft > o fx- 

f f 

By the maximum principle, Lemma 2.3 in 0, one has 

mimp < tp, for (t, x) G £Xr, (2.7) 

where tt T = [-1, 1] x [0,T] and <9'fi T = dtt T \{{x,T)\ - 1 < z < 1}. Note that / and g 
satisfy (12. ip . we have 

ip = Cu x = l + u x = l + F + g = f + g>e 1 , at x = 1, 

for some ei > 0. The oddness of the function u) yields 

<p(-l,t) = <p(l,t). 

Therefore, the estimate (12.71) is equivalent to 

up > min{l, min{/(l, i) + g(t)}} > min{l, = e 2 > 0. (2.8) 

Define TV = ^ 2 /2 and 

/i(t) = maxjmax F„(i,t),0} and V(t) — sup — + / h(s)ds. (2.9) 
ze[o,i] a'n T 2 J 
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By straightforward computation, the equation (12.61) can be rewritten as 

W t = \{W- 1/2 W X ) X - {tp x - F x )—W x - F xx . (2.10) 
It follows from ([22]) and f[2~T0D ) that 

(W - V) t < ^{w-y\W - V)^ - (cp x - F X )^(W - V) x , in Q, 
W-V<0 on &Q, 
where we used V x = 0. By the maximum principle again, we have 

sup (p = sup V2W < y/2V(T). (2.11) 

As long as ip is positive and bounded, the equation (II. 5p is uniformly parabolic. Therefore, 
the global existence of the problem (I1.5P is a consequence of Theorem 12.14 in [9]. 

Set \1/ = if — F — g. If f(x,t) = f(x) and g(x) = g for t > T , then equation (12. 6p becomes 

•"(£). (2 ' 12) 

for t > T . Note that ^(-1, t) = t) = 0, multiplying the both side of fl2TT2|) with ^ and 
integrating by parts, one has 

r 1 r 1 i 2 

d t / \^\ 2 (x,t)dx + 2 / J— ^-(x,t)<te = 0. (2.13) 

By the maximum principle, for t > T , the solution ^ of equation (I2.12p admits the 
following estimate 



* < max{ sup *(x,T o ),0}. 

-Ki<l 



This implies 



sup if < sup (F + g) + max{ sup \l/(x, T ), 0}. 

-1<z<1, -l<a<l -1<X<1 

t>T 

Combining with (12. lip gives 

sup <p < -L (2.14) 

-l<a;<l, £Ol\ 
t>T 

for some constant «i > 0. Thus, taking (12. 8p and (12.141) into account, it follows from (12 . 13[) 
that we have 



d t J \^\ 2 (x,t)dx + a! J \^ x \ 2 (x,t)dx < 0. 
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Using Poincare inequality gives 

d t j \^\ 2 (x,t)dx + 2a j \^\ 2 (x, t)dx < 
for some a > 0. This implies 



Furthermore, the Nash-Moser iteration (Theorem 6.17 and Theorem 6.30 in [9]) shows that 
that there exist to and C (independent of t) such that 



Using estimate (12. 151) yields that 

sup \^{x,t)\ < Ce~ at . 

-l<x<l 

This finishes the proof for part (a). 

(Part b). Let us start with the equation (12. 6j) . If F(x,t) = f(x) and set 

tp(x,t) = ip(x,t) - F, 

then ib satisfies the equation 




(2.15) 




= ~^XX ylPi 

if 

By the maximum principle, Lemma 2.3 in [9], one has 



min-0 < ib < max 
9'n d'n 



Thus 



min ib + F < ld < max ib + F. 

d'U d'Q 



Note that at x = 1, 



ib = ip-F= 1 + F + g-F 



i + g, 



so 



min^(l,t) + = 1 + min#(t) + /(l) 



1 = /(l) + min 5-. 
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When t = 0, 

min(^ + F) = min(l - (/ - 1)) + (/ - 1) = 1 + (/ - max/). 
If (12.41) holds, then there exists an e > such that 

max / — min f < 1 — e. 

Hence tp > e. Similarly, 

maxtp + F = max{/(l) + maxg, 1 + / — min /}. 

<9'f! 

Hence, we can also show global existence of the problem (11.51) via Theorem 12.14 in [9]. 

Similar to part (a), we can get the asymptotic behavior of the solution by energy estimate 
and Nash-Moser iteration. This finishes the proof for part (b). □ 



3. Numerical studies of existence of biological solutions 

The assumptions on / and g in Theorem 12.11 are sufficient conditions to guarantee global 
existence of solutions, thus it would be interesting to examine the solutions when these 
assumptions are violated. In this section we utilize a numerical method to find approximate 
solutions of the nonlinear problem. 



3.1. Numerical method. We consider solving a more general nonlinear problem 

\ dx ) dt dx \dx ^ dx X ' 
u(0,f) = 0, ?j± + p?ji-(f(x,t)-l)=g{t), (3 ' 1} 
k u(x,0) = 0. 

with a finite element method in the Sobolev space if 1 (0,1), that is, the square-integrable 
functions up to the first order weak derivative. 

Choose a uniform time step At > and denote the time points when solutions are sought as 
t k = kAt, k = 0, 1, • • • . The numerical scheme is to find u k+1 (x) £ H 1 ^, 1) with u k+1 (0) = 0, 
such that for any test function £ -ff 1 (0, 1) with 0(0) = 0, 

i i 



dx 



At 



du k 



dx 



d u k+1 - u k 



dx 



At 



(/ 



k+l 



(3.2) 
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The spatial domain [0, 1] is uniformly discretized into n equal-sized sub-intervals with mesh 
size h — -, and mesh points are denoted as X{ — {i — l)h,i = 1, ■ ■ • , n + 1. The space 
if 1 (0, 1) is approximated by P 1 finite element space: 

14(0, 1) = {v h e C°(0, 1) : Vh is a linear function on each subinterval [xi, i — 1, ■ • ■ , n}. 
Numerical tests show this scheme is first order accurate in time (data not shown). A Matlab 



version of the code has been provided in the website http: //www.cst. cmich.edu/users/zhenglx/ 
In all the numerical simulations in this work, we have chosen n = 1000 and At = 10~ 5 , and 
each simulation result is almost identical to that with the more refined choice n = 2000 and 
At = 5 x 10~ 6 . 

3.2. Comparisons between friction and viscosity. The analysis in Appendix lAl suggest s 
that the viscosity and friction have an additive effect on the time scale of capillary extension, 
and the viscosity of biological values is negligible compared with the friction. In this section, 
we use numerical simulations to test the effect of viscosity on solutions with respect to 
friction. 

First, we fix = 0.01, g = 5.7, / = 1, and compare the solutions in three cases: inviscid 
case, pL — 0; biological 10 4 (corresponding to dimensional viscosity \x = 10 4 — f ); 

and exaggerated case, /x = 1 (corresponding to biological viscosity jl = 10 8 comparable 
to that of pitch at 20°C). The numerical results are shown in Fig. 13.11 The fi = (results 
not shown) is almost identical to the \i = 10~ 4 case and in both cases the solutions reach the 
same steady state at roughly the same time t = 0.1. In the exaggerated case, the solution 
reaches the steady state at t — 5. Therefore, when the viscosity is significantly larger than 
friction, the solution goes to the steady state much slower. 

Second, we change the friction to /3 = 1 and re-do the above three simulations. The 
numerical results are shown in Fig. 13. 21 The solutions in the inviscid and biological cases 
are also not distinguishable, and both reach the steady state at t = 7, while the solution 
in the exaggerated case gets to the steady state at t — 15. Although the \i — 1 case has 
the decreased capillary extension speed, it is not significantly different from the inviscid case 

(Fig.E3[b]). 

These comparisons confirm the suggestion at the beginning of this section, that is, when 
we focus on the biological values, it is valid to neglect the viscosity. 



3.3. Examples examining Theorem l2.ll (a). In all the examples of this subsection, we 
let j3 = 1, n = 0, and g = 0. 
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Figure 3.1. Solutions when (3 = 0.01, g = 5.7,/ = 1. [a]: biological case, 
\x = 10 -4 . The solution curves are at time =0.01, 0.02, ■ • ■ ,0.1 (bottom to top), 
[b]: exaggerated case, /i = 1. The solution curves are at time=0.5, 1, ■ ■ • ,5 
(bottom to top). 



6r 6r 




Figure 3.2. Solutions when (3 = I, g = 5.7, and / = 1. In both [a] and [b], 
the solution curves are at t — 1, 2, ■ ■ ■ , 7 (bottom to top), [a]: biological case, 
H = 10~ 4 . [b]: exaggerated case, /i — 1. 



First, we choose / = — 10a; 2 + 10 + 10~ 6 . This choice satisfies the conditions in Theorem 12.11 
(a). The numerical solution from t — 1 to t — 10 is illustrated in Fig. l3.3[ a]. It is clear that 
1 + u x is non-negative for all time, and the solution converges to u s = — yi 3 + (9 + 10~ 6 )x, 
the steady state solution. The maximum error between the numerical solution and u s is 
shown in Fig. l3.3| b]. and it tends to zero in time exponentially. 

Second, we choose / = 10x 2 + 10. The choice of / violates the condition f xx < in 
Theorem l2.ll (a). The value of 1 + u x is negative for x G [0, 0.015] when t = 0.03 (Fig. l3.4j) . 
A few time steps later the numerical solution blows up near x = 0. 
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Figure 3.3. [a] Solution curves (solid) and (1 + u x ) curves (dashed) in the 
case of p = 1, y, = 0, g = 0, / = -10x 2 + 10 + 10" 6 , at time t = 1, 2, ■ ■ • , 10 
(both from lower to upper), [b] \v,h — u s \ 00 with respect to time. The y-axis is 
in log scale. 



0.2 

0.15 

= 0.1 

oi 0.05 
ui 

0J 
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| -0.05 
3 -0.1 
-0.15 
-0.2, 







0.005 0.01 
x, Lagrangian coordinate 



_1+u 



0.015 



Figure 3.4. Numerical solution and 1 + u x at t = 0.026, 0.027, • • • , 0.03 (up- 
per to lower) in the case of /3 = 1, \i — 0, g — 0, / = 10x 2 + 10. 



Notice in the mathematical model the density serves as a pressure exerting in the direction 
— V/, that is, from high pressure regions to low pressure regions. The breakdown of the 
second example stems from the high density at the tip driving towards the root. At the root, 
the boundary condition u = makes the root cell nowhere to escape. So even after it is 
smashed to zero length, the pressure still cannot be balanced, which results in the breakdown 
of the model. In the first example, the pressure with the same magnitude of gradient exerts 
on the system but towards the tip direction. The free boundary condition at the tip releases 
the pressure and avoid the breakdown of the system. 

These two examples implies that, if the pressure is monotonic along the vessel, the position 
of high pressure relative to the fixed boundary is important for the existence of biological 
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solutions. But positioning the high pressure at the tip does not necessarily break down the 
system. For instance, we also tested / = 1 + x 2 , in which case the high pressure is at the 
tip but the global biological solution still exists. Comparing with the case of / = 10 + 10x 2 , 
which is of higher pressure drop from the tip to the root, we deduce that the pressure drop 
from tip to root must be large enough to break down the system. 

3.4. Examples examining Theorem 12.11 (b). The same as the previous section, we also 
set (3 = 1, /j, = 0, and g = for all the examples in this subsection. 

First, let / = 0.4999 cos(100x) + 0.5. This choice satisfies the conditions in Theorem l2.ll 
(b). The numerical solution from t — 0.1 to t — 1 (Fig. 13 .5 [ a]) indicates that (1 + u x ) is 
non-negative for all time, and the solution converges to 0.004999 sin(100x) — 0.5a;, the steady 
state solution. The maximum error between the numerical solution and the steady state 
solution is shown in Fig. l3.5[ b], and it tends to zero in time exponentially. 




Figure 3.5. [a] Solution curves (solid) and (1 + u x ) curves (dashed) in the 
case of p = 1, n = 0, g = 0, / = 0.4999 cos(100a;) +0.5, at time t = .1, .2, • • • , 1 
(from upper to lower), [b] \uh — u s \ 00 with respect to time. The y-axis is in 
log scale. 



Second, we change to / = 2 cos (24a;) +2.1. This choice violates the condition max^^ii f{x) — 
min xg [ 0> i] f(x) < 1 in Theorem 12.11 (b). The numerical solution at t — 10~ 3 is illustrated in 
Fig. l3.6t which shows that 1 + u x is negative at several points, where the solution blows up 
later. 

Note that the pressure gradient in the blow-up case is even smaller than that in the non- 
blow-up case, therefore, the blow-up in the second example cannot be attributed to the high 
pressure gradient. However, the pressure amplitude in the blow-up case is larger than the 
non-blow-up case. Therefore, it is pressure's amplitude that breaks down the system. 
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Figure 3.6. Solution curve (solid) and (1 + u x ) curve (dashed) in the case 
of p = 1, n = 0, g = 0, / = 2 cos(24x) + 2.1, at time t = 10" 3 . 

4. Numerical simulations of angiogenesis 

In this section, we apply the mathematical model to simulate angiogenesis under three real 
biological conditions: capillary extension without proliferation, extension with proliferation, 
and capillary regression. 

4.1. Capillary extension without proliferation, but with variable friction. In the 

experiments of [T4] . the endothelial cells are given X-ray irradiation. With enough doses of 
irradiation, DNA synthesis are stalled and cells lose the ability to proliferate. This means the 
cell density will keep the same from the beginning. Therefore, the cell density is always equal 
to the reference density, that is, f(x, t) = 1 for all the points and all the time (see Appendix 
for non-dimensionalization). The viscosity and protrusion force are chosen as \x = 10~ 4 and 
g = 4.7, within the biological ranges as discussed in the Appendix. However, the friction is a 
variable in time. Because, at the initial stage of angiogenesis, endothelial cells have to over- 
come the tight association with surrounding cells and basement membrane. These obstacles 
are rapidly reduced after the initial stage, so the friction is smaller with the progression of 
angiogenesis. Therefore, we assume the friction is a time-dependent function 



whose graph is plotted in Fig. l4.1[ a]. This is a pretty rough approximation of the true 
variation of the friction. Indeed the friction depend on many other conditions, such as the 
proteins regulating EC-environment adhesion, association with pericytes, etc. Readers are 
referred to [6] for more details. 



P(t) = 0.01 + lOOe 



-1.8* 
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The numerical solutions are shown in Fig. l4.1[ b]. At time t — 1, only the points near 
the tip have observable motion, because the friction is very large initially and only the cells 
near the tip can be dragged with remarkable distance by the protrusion force exerted at the 
tip. Up to t = 3, all points migrate toward the tip direction. The deformation of the tip 
is 3.9 at t — 4 and 4.7 at t — 7. These correspond to the whole capillary length 490 fim at 
day 4 and 570 fim at day 7 after the initiation of angiogenesis, very close to the rat corneal 
experimental results in [T4] . 




FIGURE 4.1. Extension of capillary without proliferation, f(x,t) = 1. [a] 
Friction as a function of time, (3 = 0.01 + 100e -1 ' 8 *. Y-axis is in log scale, [b] 
Deformation curves at time t = 1,2, ••• ,7 (lower to upper). The curves of 
t — 5, 6, 7 almost coincide. Other parameters: g = 4.7, /x = 10~ 4 . 



Although there is no cell proliferation, the protrusion force of the tip cell drags the capillary 
to extend. The ultimate deformation is limited by the elasticity and is given by the steady 
state solution u s = g/E. For the above example, this value is 4.7, which has been achieved 
by the numerical simulation. However, the lack of proliferation does affect the capillary ex- 
tension: the capillary only reaches 570 /zra, but the distance from the root to the chemotactic 
source is about 2 mm in the experiments of [14] . The limited capillary extension without 
proliferation is an important feature of biological growth, but it is very hard to capture using 
reaction-diffusion-convection models. For example, many reaction-convection-diffusion mod- 
els (e.g. [TO] ) only tracks the motion of the tip cell, which is modeled as migrating towards 
the chemotactic source without any restriction of trailing cells. Thus, it falsely predicts the 
capillary reaches the chemotactic source even when there is no cell proliferation. 

4.2. Capillary extension with proliferation and variable friction. Comparing with 
the last example, we change the cell density to / = 2t + 1. The biological meaning of this 
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function is that the cells uniformly proliferate in the whole capillary: each cell generates 
two more cells per day. The numerical results are shown in Fig. 14.21 The value of u of the 
capillary tip is equal to 0.53 at t — 1, 8.64 at t — 4, and 18.5 at t — 7. These represents the 
capillary reaches 63 fim at Day 1, 964 /im at Day 4, and 1950 fim at Day 7. These results 
reproduce the data in the experiments of [2] and at Day 7 the capillary has already reached 
the chemotactic source. 




x, Lagrangian coordinate 

Figure 4.2. Capillary extension with cell density / = 2t + 1. Other param- 
eters: j3 = 0.01 + lOOe -1,8 *, g = 4.7, n = 10~ 4 . Curves from lower to upper: 
t = 1, 2, • • • , 7 days. 

In this example, both the tip cell protrusion and stalk cell proliferation contribute to the 
capillary extension. Comparing with the last example, cell proliferation is more decisive in 
supporting the extension of the capillary until it reaches the chemotactic source. Mathemat- 
ically, the increase of cell density alters the balance between the previous stress and density. 
The cells enlarge the deformation to adjust to the higher density, which allows the extension 
of the vessel. Therefore, the tip cell migration plays a leading role in capillary extension 
while stalk cell proliferation provides necessary material to support the extension. 

4.3. Capillary regression. The newly formed blood vessels in pathological corneal angio- 
genesis can be treated by the drug Bevacizumab as in [TTJ to induce apoptosis of endothelial 
cells, thus the regression of blood vessels. To simulate this therapy, we reset the friction as a 
constant /3 = 1, and the cell density as / = max(l — tx, 0) (see Fig l4.3| a]). which represents 
the cells begin to die from the tip to the root. Because the tip cell is dying, it loses the ability 
to produce the protrusion force. Therefore, g = 0. The numerical results in Fig J4.3[ b] show 
that the deformation of the capillary tip is —0.9 at t = 5, which means the capillary shrinks 
by 90 \xm after 5 days' regression, so only 10 fim is left. The deformation gradient 1 + u x 
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at t = 5 is zero for x G [0.2, 1], indicating the cells initially in this interval have died out. 
This simulation basically captures the experiments observed in human corneal angiogenesis 
treatment results in 




Figure 4.3. Capillary regression, [a] Cell density / along the capillary at 
time t = 1,2,3,4,5 (from upper to lower), [b] Deformation u curves (solid) 
and 1 + u x curves (dashed) are at t — 1, 2, 3, 4, 5, both from upper to lower. 



Contrary to the last example, the cell density decreases with respect to time in capillary 
regression. In this model, the stress/stain is supported by the density. Therefore, the loss 
of density induces the decrease of stain, thus the regression. This phenomenon is almost 
impossible to be modeled by reaction-convection-diffusion models because they lack the 
mechanisms to retract the capillary. For example, in the model of [10] where only the tip 
cell is modeled, the tip cell would stay where it used to be in space once the chemotactic 
cues are removed by therapy, thus regression observed in experiments would not occur. 

5. Conclusion 

This article analyzes a one- dimensional blood vessel growth model with both analytical 
and numerical tools, and demonstrates its power in simulating various biological growth 
cases. 

Theorem l2.ll concludes that the global biological solution exist if (a) the cell mass density 
is a decreasing function along the capillary from the root to the tip, or (b) the amplitude 
of the cell density along the capillary is less than the unity. The system may break down if 
these two conditions are violated. Indeed numerical examples indicate that if there is a big 
pressure drop from the tip to the root direction or a big oscillating amplitude of density, then 
the solution would probably blow up at the low density point. The underlying mechanism 
of breakdown is the failure of linear elasticity to resist the high pressure impact. 



VISCOELASTIC ANGIOGENESIS 17 

Besides modeling the biological growth, this current model can also be used to describe the 
slow deformation of one-dimensional thermoelastic bar, where the stress tensor is the same 
as (11.21) except the cell density is interpreted by temperature (cf. [12] )■ Therefore the analysis 
can be applied in the thermoelastic case. For example, if the temperature decreases from the 
root to the tip of the bar, or has small oscillations along it, then the global solution exists. 
However, if the tip of the bar is overheated or a large temperature oscillation is applied on 
the bar, then the bar may be broken by the expansion pressure produced by temperature 
differentials. 

The reaction-convection-diffusion models are predominant in continuous models of angio- 
genesis, but they encounter great difficulties in simulating capillary growth without prolif- 
eration and capillary regression. The difficulties come from the fact that they are unable to 
build connections between tip cell protrusion, change of stalk cell density, alteration of the 
stress, and capillary extension. In contrast, the model of [6] regards the whole capillary as 
one viscoelastic material dragged by the tip cell and directly relates cell density with elastic 
deformation. Plus the nonlinear interactions with surroundings, this model is capable of 
simulating various capillary growth patterns. 

Appendix A. Mathematical model, non-dimensionalization, and parameters 
The blood vessel capillary extension is described in the reference coordinate x G [0, L] by 

2P\l + u x \ . da 

u = 7-, (A.l) 

r ox 

a = £* +A ^_fc*>, (A.2) 
ox ox po 

u(0,t) = 0, a(L,t)=G, (A.3) 
u(x,0) = 0, (A.4) 

(A.5) 

where L is the initial capillary length, u(x, t) is the displacement of material point x at time 
t, (3 is the friction, r is the capillary radius, E is the Young's modulus of endothelial cells, 
\x is the dynamic viscosity, po is the initial cell density, F is the cell density multiplied with 
1 1 + u x | , and G is the protrusion force at the tip of the capillary. For the derivation of this 
model and other related models, see [6]. 

Choose L as the characteristic length, which is 100 /im based on the experiments in |14j . 
Because angiogenesis typically runs a couple of days from the initiation to the penetration 
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of tumor/tissue, we choose T = 1 day as the characteristic time scale. Let x = x'L,u = 
u'L, t = t'T, and insert these relations to the above system, where we then delete the primes 
for simplicity. Therefore we obtain for x G [0, 1], 

B\l + u x \u = (A.6) 

o = + (A.7) 

OX ox 

u(0,t) = 0, a(l,t) = g, (A.8) 

u(x,0) = 0, (A.9) 

where /3 = ^ = aiid «, = §. 

The Young's modulus E for endothelial cells is between 1.5 ~ 5.6 x 10 3 according; to 
[TJ. The viscosity /i is not available for endothelial cells, so we replace it with the value for 
fibroblasts [15] : a = 10 4 ^ (comparing with 10" 3 of water, 3 x 10 4 ^ of tar, and 
2.30 x 10 8 ^| of pitch , all at 20°C [16]). The estimate of (3 will be highly dependent on 
the cell-environment contacts. For example, in the case of bovine aortic endothelial cells 
(BAECs) spreading on polyacrylamide gels, the friction is deduced to be /3 = 10 3 [5]. 
However, ECs experience much stronger friction or resistance in the in vivo condition, because 
ECs are in association with surrounding cells such as pericytes. Therefore, we choose a large 
range for the friciton: 2.5 x 10 3 ~ 6.5 x 10 6 ■ The radius of blood vessel capillary is about 
10 jjm. The protrusion force G is about 10 4 -j^ as measured in [13]. With these values, the 
non-dimensionalized parameters become \i = 10 -4 , (3 G [0.01, 100], and g G [1.7, 6.7]. 

Both the viscosity and friction resist the motion of capillary, and their relationship can be 
analyzed from a model problem. Assume the deformation u(x,t) is a periodic function in 
the free space with period 2ir, satisfying u(0,t) = u(2ir,t) = 0, and 

Pik = u xx + fiu txx . (A. 10) 

oo 

Express u(x, t) in the Fourier series Uk(t) sin(for), and insert it to flA.lOj) . Then we obtain 

fe=i 

(P + [ik 2 )u Kt = -k 2 u k (A.ll) 



for each mode number k > 1. The time scale for each mode to reach the steady state is 
P/k 2 + fi. Therefore, the time scale for the whole system is P + \i. Given the above values 
of \i and P, the viscosity is far less than the friction. This implies the viscosity would have 
negligible effects compared with the friction on the capillary extension. 
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